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Time-accurate, Euler, quasi-ID and 2D-axisymmetric computational fluid dynamic 
(CFD) simulations of a three-stream plug nozzle were performed as part of an effort to 
provide an efficient and accurate CFD-based approach to modeling nozzle dynamics. The 
CFD code used for the simulations is based on the space-time Conservation Element and 
Solution Element (CESE) method. Steady-state results were validated using the Wind-US 
code and a code utilizing the MacCormack method while the unsteady results were partially 
validated via an aeroacoustic benchmark problem. The CESE steady-state flow field 
solutions showed excellent agreement with solutions derived from the other methods and 
codes while preliminary unsteady results for the three-stream plug nozzle are also shown. 
Additionally, a study was performed to explore the sensitivity of gross thrust computations 
to the control surface definition. The results showed that most of the sensitivity in 
computing the gross thrust is attributed to the control surface stencil resolution and choice 
of stencil end points and not to the control surface definition itself. Finally, comparisons 
between the quasi-ID and 2D-axisymetric solutions were performed in order to gain insight 
on whether a quasi-ID solution can capture the steady and unsteady nozzle phenomena 
without the cost of a 2D-axisymmetric simulation. Initial results show that while the quasi- 
ID solutions are similar to the 2D-axisymmetric solutions, the inability of a quasi-ID 
simulation to predict two dimensional phenomena limits its accuracy. 
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I. Introduction 

A s research and development budgets become more constrained with time, wind-tunnel model testing becomes 
less affordable for aircraft design. Computer-based models have the potential to enable rapid aircraft design 
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and analysis at reduced expense compared to wind-tunnel model testing. This is exceptionally true in the early and 
intermediate stages of the design phase. As such, there are efforts in the aerospace industry to develop models that 
maximize accuracy while minimizing the turnaround time and the computational resources needed. The model can 
encompass the entire aircraft, a particular system, or even a subcomponent of a system. This paper presents analyses 
of a simulated three-stream plug nozzle and investigates both steady and unsteady nozzle flow behaviors. The 
analyses encompassed quasi- ID and 2D-axisymmetric time-accurate Euler simulations. One of the primary goals of 
the simulations was to develop an optimized computational fluid dynamic (CFD) -based nozzle model that could be 
used as part of a system level model. In this paper, optimization is defined as minimizing the run time and the 
amount of computational resources needed to compute the time-varying gross thrust from a nozzle while retaining 
an accurate prediction. This paper compares the results of the simulations and also explores factors involved in 
calculating gross thrust for a CFD -based nozzle model. 


II. Geometry and Numerical Modeling 

The nozzle chosen for the CFD simulations consists of three internal ducts with a downstream conical plug or 
center body. For nomenclature clarity the ducts are denoted as follows (in increasing radial order): inner duct, 
primary duct, and noise shielding (or simply shield) duct. The simulations used quasi- ID and 2D-axisymmetric 
grids to represent the nozzle. Due to the proprietary nature of the actual nozzle, this report will only show a generic 
layout of the nozzle, shown in Fig. 1, and present results in a non-dimensional format based on the freestream 
density, speed of sound, and a characteristic length. 



Figure 1. Generic three-stream nozzle geometry (axisymmetric cut). 

The space-time Conservation Element and Solution Element (CESE) 1,2,3 method, originally proposed by Chang 
at NASA’s Glenn Research Center, was the primary numerical method chosen for the CFD simulations. The CESE 
method is a time accurate formulation with flux-conservation in space-time and utilizes a staggered mesh. The 
staggered mesh allows for the scalar quantities (e.g. pressure, density) to be computed at the cell center while the 
vector flow quantities (e.g. momentum, velocity) are computed at the cell faces. One benefit of this is that it 
eliminates the need to solve the Riemann type problem along the cell interfaces. Further, the CESE method is 2 nd - 
order accurate in both space and time with the stability condition of Courant-Friedrichs-Lewy (CFL) number less 
than one. In the past, the CESE method has been applied to a wide variety of CFD and computational aeroacoustic 
(CAA) problems 4,5,6 . The CESE method has been proven to be a robust and efficient tool for dynamic simulations 
and aero-acoustic computations. 

In addition to the CESE method, a code utilizing the MacCormack scheme 7,8 and the Wind-US code 9 were 
chosen to help validate the CESE derived steady-state flow fields. The MacCormack scheme was chosen for its 
simplicity in modeling a quasi- ID flow field. However, the chosen MacCormack scheme -based code only allows 
for quasi- ID simulations. The Wind-US code was then chosen for comparing the 2D -axisymmetric simulations and 
has been shown to model nozzle flows effectively 10 . Wind-US is a general purpose compressible flow solver 
developed by the NPARC Alliance that can run in either Reynolds -Averaged Navier-Stokes (RANS) or Euler 
modes. The default scheme for Wind-US is Roe’s 2 nd -order up-winding scheme 11,12 and unless otherwise specified, 
was run in the Euler mode. 


III. Flow Field Correlation and Validation 
A. Steady-State Quasi-ID Simulations 

The CESE quasi- ID Euler code was used to compute a steady-state solution for each duct of the nozzle. Each 
duct was represented by its own 200 uniform- spaced element grid. Two additional solutions were computed for 
comparison and validation: a quasi- ID, perfect gas, isentropic analytical solution and an independent quasi- ID 
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solution utilizing the MacCormack scheme, both provided by Mr. Joseph Connolly. Results from these simulations 
for the primary duct are shown in Fig. 2. It can be seen that the CESE solution is in excellent agreement with the 
MacCormack and analytical solutions for the steady-state flow field. 




Figure 2. Quasi-ID pressure (left) and u velocity (right) axial distributions for the primary duct. 


B. Steady-State 2D-Axisymmetric Simulations 

2D-axisymetric simulations were performed to assess the limitations of the lower resolution quasi- ID 
simulations. Two codes were used: the 2D CESE Euler code and Wind-US. The CESE simulations used a constant 
time step while the Wind-US simulations used a constant CFL number. For the Wind-US simulation, a 129,000 cell 
structured grid modeled the external plume and plug regions as well as all three internal ducts of the nozzle. The 
CESE simulation used a diagonalized derivative of the Wind-US grid, resulting in a 258,000 cell triangular grid. An 
example of the diagonalization process is shown in Fig. 3. 



Figure 3. Example of diagonalizing a structured grid. 


A grid resolution study was performed on three levels of refinement of the Wind-US grid. The coarse level used 
every fourth grid point in each computational direction, the medium level used every other grid point in each 
computational direction, and the fine level used all grid points. The mass flow rate through the mixing region of the 
inner and primary ducts was computed via the Wind-US utility code CFPOST at each grid level, the results of which 
are shown in Table 1. Note that the percent difference between the coarse and medium grid levels (c-m) was 
defined as: 


% Difference = 100 x 


iqq — iqq 

coarse medium 
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The percent difference between the medium and fine grid levels (m-f) was defined as: 
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while the percent mass flow gained was defined as: 

iqq — iqq 

r\ / 7i jf /— r j < r\r\ Zone6,i=last Zon€5-\-Zone6 ^ _ x 

%MassFlow = 1 00 x — (3) 

iqq 

n Zone5+ZondS 

Figure 4 shows a graphical representation of the grid zones presented in Table 1. It can be seen from Table 1 that 
the fine grid level is reaching grid convergence and is sufficiently converged for the presented analyses. Assured of 
grid convergence, u velocity and pressure through the primary duct are shown in Fig. 5 for the 2D-axisymmetric 
Euler simulations. Note that the ID cuts shown in Fig. 5 were taken as the streamline originating at the duct inlet, 
center height, to reduce the dimensionality for later comparisons with the quasi- ID solution. The streamline 
approach was used over an averaging technique to better capture the presence of shock waves. From Fig. 5 it is 
shown that the CESE solution agrees well with the Wind-US solution (within 0.407% for pressure and 0.163% for u 
velocity). Both solutions show a sharp increase in pressure (with a corresponding decrease in u velocity) near x=0.9, 
which is due to an oblique shock wave off one of the duct splitters. However, the CESE solution predicts the shock 
location just upstream of the Wind-US solution, resulting in a higher minimum pressure and a lower maximum u 
velocity. 


Table 1. Mass flow rates at various grid levels for the Euler Wind-US grid. 


Zone 

Coarse 

Medium 

Fine 

% Difference (C-M) 

% Difference (M-F) 

Zone 5 + Zone 4 (Mast) 

2.413 

2.401 

2.398 

0.500 

0.125 

Zone 6 (i=l) 

2.379 

2.393 

2.401 

-0.585 

-0.333 

Zone 6 (Mast) 

2.400 

2.401 

2.406 

-0.042 

-0.208 

% Mass Flow Gained 

-0.539 

0.000 

0.334 

- 

- 



Figure 4. Graphical representation of the grid zones used in Table 1. Note that the number is the zone 
number. 



Figure 5. 2D-axisymmetric pressure (left) and u velocity (right) axial distributions for the primary duct. 
C. Dynamic Quasi-ID Simulation 
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For ease of validating an Euler, quasi- ID, unsteady flow field, the CESE method was applied to an aeroacoustic 
benchmark problem 13 rather than directly to the nozzle of interest. The aeroacoustic benchmark problem simulated 
a downstream running acoustic wave at the inlet of a convergent-divergent nozzle with the acoustic wave 
disturbance defined as: 
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Note that the benchmark problem was poised in a non-dimensional format and the mean flow boundary conditions 
are shown in Table 2, per reference 13. Although the problem is over specified for a subsonic inflow, it is still a 
useful benchmark problem with steady and unsteady solutions shown in Fig. 6 and Fig. 7. 

Table 2. Non-dimensional mean flow bo undary conditions fo r the aeroacoustic benchmark problem 


M Met 

0.2006533 

Pinlet 

1/y 

Pinlet 

1 

Pexit 

0.6071752 


CESE (20D) 



Figure 6. Steady-state pressure distribution. Note that the numbers in parenthesis are the number of grid 
elements. 
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CESE (200) CESE (200) 



Figure 7. Dynamic pressure distribution (left) and time history at the nozzle exit for one cycle (right). 

Note that the numbers in parenthesis are the number of grid elements. 

The results, computed using both 200 and 400 uniformly-spaced element meshes, show excellent agreement with the 
analytical solutions provided for the benchmark problem and serves as confirmation of the results shown by Wang 
et. al. 14 , who previously demonstrated success in modeling the quasi- ID steady-state flow field as well as the time- 
varying dynamics of the aeroacoustic benchmark problem with the CESE method. While this validates the CESE 
quasi- ID code, the authors are not aware of an equivalent benchmark problem to analytically validate a 2D- 
axisymmetric unsteady simulation. 


IV. Comparisons Between Quasi-ID and 2D-Axisymmetric Simulations 
A. Steady-State Flow Field 

After validating the CESE steady-state solutions, the 2D-axisymmetric solution was compared to the quasi- ID 
solution, shown in Fig. 8, to identify any differences in the computed flow fields. The quasi- ID solution agrees well 
with the 2D-axisymmetric solution up to about x=0.7, with the quasi- ID solution slightly over predicting the 
upstream u velocity (and subsequently under predicting the upstream pressure). However, the quasi- ID solution 
fails to pick up the oblique shock wave computed in the 2D-axisymmetric solution, which is a primary deficiency of 
a quasi- ID simulation. The oblique shock wave is represented by the sharp increase in pressure and decrease in u 
velocity at x=0.9 and is the main reason for the discrepancy between the solutions downstream of x=0.7. 



Figure 8. Pressure (left) and u velocity (right) axial distributions for the primary duct. 
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B. Gross Thrust Computation 

While examining flow field quantities such as pressure and u velocity can provide insight into the differences 
between the quasi- ID and 2D-axisymmetric solutions, further insight can gained by examining the computed gross 
thrust for the nozzle. Gross thrust is of particular interest as it is one of the main parameters used when 
designing/sizing a nozzle for integration with an aircraft propulsion system. Analytically, the gross thrust is defined 
over a control volume as: 


F = J pu 2 dA + J (p - p x )dA 


( 8 ) 


However, the gross thrust can be numerically computed by defining a control surface with a discrete stencil and 
numerically integrating the thrust along the stencil. For 2D-axisymmetric geometries, Eq. 8 becomes: 
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One potential flaw in Eq. 9 is that it will only compute quantities over the discretized area, which depending on how 
the stencil endpoints are chosen, might not represent the true area. This can be caused by having the stencil start and 
end nodes not on the true geometrical boundary. To mitigate the effect of the potential gaps at the nozzle walls, the 
computed gross thrust is corrected by the known geometrical area. 


F = 


A 

, actual 


A 


stencil 


(ii) 


For computing the gross thrust, three different control surfaces were used and are shown in Fig. 9 through Fig. 11. 
Because Eq. 9 and Eq. 10 are specifically for 2D-axisymmetric geometries, Eq. 8 was computed directly at the 
nozzle exit for the quasi- ID solutions, with a representative control surface shown in Fig. 12. 



Figure 9. 2D-axisymmetric gross thrust control surface (shown in red), Configuration 1. 



Figure 10. 2D-axisymmetric gross thrust control surface (shown in red), Configuration 1.5. 
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Figure 11. 2D-axisymmetric gross thrust control surface (shown in red), Configuration 2. 



Figure 12. Quasi-ID gross thrust control surface (shown in red). 

After developing the means to compute the gross thrust with different control surface definitions, a sensitivity 
study to the stencil resolution was performed. As a reminder, the stencil is the discrete points that represent the 
control surface, with a generic stencil concept shown in Fig. 13. 



Figure 13. Sample control surface stencil. Note that the orange squares represent the discrete points along 
the control surface (shown in red). 

Three separate 2D-axisymmetric steady-state solutions were used for the stencil resolution study: the Euler CESE 
and Wind-US solutions presented earlier and a 2D-axisymmetric viscous Wind-US solution which used the 
Menter’s Shear Stress Transport (SST) 15 turbulence model. The viscous Wind-US solution, denoted as Wind-US 
SST, was simulated to gain an insight to see if adding viscosity would reduce some of the convergence and 
conservation errors. The Wind-US SST simulation used a modified grid comprised of 401,000 structured cells with 
a max y+ of 0.49. As a check, a grid resolution study was performed on the viscous grid in the same fashion as the 
Euler grid, the results of which are shown in Table 3. Just like the Euler grid, the viscous grid is shown to be 
approaching grid convergence and shows sufficient grid convergence at the fine level for the studies presented. 


Table 3. Mass flow rates at various grid levels for the viscous Wind-US grid. 


Zone 

Coarse 

Medium 

Fine 

% Difference (C-M) 

% Difference (M-F) 

Zone 5 + Zone 4 (i=last) 

2.412 

2.391 

2.384 

0.878 

0.294 

Zone 6 (i=l) 

2.412 

2.391 

2.384 

0.878 

0.294 

Zone 6 (Mast) 

2.416 

2.392 

2.385 

1.003 

0.294 

% Mass Flow Gained 

0.166 

0.042 

0.055 

- 

- 


Figure 14 shows the sensitivity to the stencil resolution for the steady-state computed gross thrust (entire nozzle) 
and the combined mass flow rate for the inner and primary ducts. It can be seen that while all the mass flow rates 
converge to within 0.28% of each other, the gross thrusts using the viscous solution converge to a lower value than 
the Euler equivalents. This is due to low momentum flow developing in the boundary-layers of the viscous solution, 
which exploits the error derived from the gap between the starting/ending nodes of the stencil and the true geometric 
boundaries. Although Eq. 10 attempts to correct for this error, it assumes that the average values over the computed 
area are roughly equivalent to the average values over the gap areas. While this is true for Euler simulations, it is 
not necessarily true for viscous simulations. Examining Fig. 14 again, it is shown that all configurations are 
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relatively converged with 60-80 nodes per stencil for the inner and primary ducts. This stencil resolution roughly 
equates to half the radial grid resolution. Although it is not clear whether this level of stencil resolution can be used 
for other nozzles flows, it will be used for the remainder of this paper. 



Combined Inner and Primary Stencil Nodes 



Combined Inner and Primary Stencil Nodes 


Figure 14. Gross thrust (left) and mass flow rate (right) sensitivities to stencil resolution. 


C. Steady-State Gross Thrust 

Table 4 shows the computed gross thrusts for various CESE Euler solutions and stencil configurations. Note that 
in the table, the percent difference was defined as: 


F — F 

%Difference = 100 x ( 12 ) 

^quasir-lD 

It can be seen that the gross thrust computed using the 2D-axisymmetric solution is within 1 .46% of the gross thrust 
computed using the quasi- ID solution. However, the gross thrusts computed using the 2D-axisymmetric solution 
are consistently lower due to 2D effects such as oblique shock waves while all are within 0.36% of each other. This 
suggests that the computed gross thrust is not overly sensitive to the presented control surface definitions. However, 
this might not be the case when computing the unsteady gross thrust and will be explored further in the next section. 

Tabl e 4. CESE-Based Computed Gross Thrust (Euler soluti ons). 


Case 

Gross Thrust 

% Difference 

Quasi-ID 

8.383 

0.00 

2D-Axi (Config. 1) 

8.266 

-1.39 

2D-Axi (Config. 1.5) 

8.290 

-1.10 

2D-Axi (Config. 2) 

8.261 

-1.46 


D. Unsteady Gross Thrust 

For the unsteady simulations, a sinusoidal disturbance was introduced at the inlet of each duct. Since the CESE 
quasi- ID solution matched the analytical solution for the aeroacoustic benchmark, a similar incoming disturbance 
and boundary condition was used. 
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where once again 
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Each simulation was run at a different discrete frequency with a constant amplitude of 0.001. The discrete 
frequencies are summarized in Table 5. Although this boundary condition had only been validated for quasi- ID 
simulations, it was applied both to the quasi- ID and 2D-axisymmetric simulations for consistency. Data from these 
simulations were then used to construct bode plots based on the input sinusoid and the unsteady component of the 
gross thrust. 


input = 0.001sm(&tf) (15) 

output =F' = F- F steady (16) 

, , . , , f output ' 

Magnitude = 20 log (17) 

^ input j 

Phase = 360/ x [t input - t outp J (18) 


Although bode plots are presented, the reader is cautioned to view them in a qualitative light rather than a 
quantitative one as there is still discussion on how to properly model a disturbance boundary condition that results in 
meaningful dynamic data. Considering this caveat, bode plots comparing the quasi- ID and 2D-axisymmetric 
simulations for the primary duct are shown in Fig. 15 (with and without the natural time delay included in the phase 
shift). The natural time delay was defined as the time it would take for a signal to travel from the primary duct inlet 
to the primary duct exit control surface. 


T = 


^ outlet ^ inlet 
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U + 

v 



(19) 


From Fig. 15 it can be seen that the 2D-axisymmetric unsteady computed gross thrust agrees well with the quasi- 
ID unsteady computed gross thrust at the lower frequencies (with the exception of the lowest frequency) but 
diverges at the higher frequencies. This is particularly true when using control surface configurations 1.5 and 2 for 
the 2D-axisymmetric computation. Recall that the CFD solution is the same for all the 2D-axisymmetric gross 
thrusts with only the post-processing stencil changing. As a check, unsteady stencil resolution studies for 
configuration 1 and 2 are shown in Fig. 16. Note that the finest stencil resolutions are those that are shown in Fig. 
15. It can be seen that unlike configuration 2, the configuration 1 control surface is still highly sensitive to the 
stencil resolution. Although this seems contradictory to what Fig. 14 shows for steady-state, it is evident that a 
significant gap exists between the stencil for configuration 1 and the true geometry which the correction described in 
Eq. 1 1 does not appear to account for in the unsteady computation. 
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Table 5. Pulsing Frequencies. 


Non-Dimensional Frequency 

0.004 

0.04 

0.1 

0.2 

0.3 

0.4 

0.6 

0.8 

1.0 

1.5 

2.0 

2.5 

3.0 

4.0 





without natural time delay (right). 






Figure 16. 2D-axisymmetric primary duct bode plot stencil resolution for Configuration 1 (left) and 
Configuration 2 (right). Note that the number of nodes is for both the inner and primary ducts combined. 
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V. Conclusions/Future Work 

Quasi- ID and 2D-axisymmetric numerical simulations were performed on a three-stream plug nozzle using the 
CESE method to aid in the development of a CFD-based nozzle flow model. First, steady-state flow field solutions 
were validated and showed excellent agreement with solutions derived from other codes and methods, 
demonstrating the capability of the CESE method. Second, a study was performed to explore the sensitivity of the 
control surface definition when computing gross thrust. Based on the study, it was shown that there was little 
sensitivity to the control surface definition itself. However, the study showed that the computed gross thrust was 
highly sensitive to the stencil resolution as well as the choice of the stencil end points. These were especially true 
for unsteady thrust computations partially due to the ineffectiveness of the stencil area correction on the unsteady 
computation. Third, the quasi- ID and 2D-axisymmetric solutions were compared. The steady-state solutions were 
shown to be similar, however, the quasi- ID solution failed to pick up the oblique shock wave at the downstream 
portion of the nozzle due to its limitations as a quasi- ID simulation. The solutions also varied when looking at the 
unsteady computed gross thrust, suggesting that the quasi- ID solution might not be sufficient for some system-level 
models. Because the disturbance boundary condition presented is most likely not suitable for real nozzle flows, 
potential future work includes the exploration of how to properly model a disturbance boundary condition such that 
meaningful dynamic data can be obtained from a CFD simulation. Also, additional unsteady simulations using the 
Wind-US code will be run to compare to the CESE-based unsteady solutions with a potentially updated disturbance 
boundary condition. 
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